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£~~je\ Method for Tracking Motion Phase of an Object for Correcting Organ Motion Artifacts in 
J\ X-ray CT Systems 



^Field of the Invention 



[001] This invention relates to medical imaging and in particular to signal processing 
methods correcting organ motion artifacts in cardiac imaging. 

Background of the Invention 

[002] X-ray Computed Tomography (CT) is a medical imaging tool providing cross-sectional 
images or images of "slices" of a human body. A patient is placed between an X-ray source and 
an array of detectors for measuring the amount of radiation passing through the patient's body. 
During the data acquisition process, the source and the detectors are rotated around the patient 
acquiring a large number of X-ray projection data. The X-ray projection data are then processed 
using a tomographic image reconstruction process in order to produce an image of a cross- 
section of the patient's body. For the data acquisition and the image reconstruction process it is 
assumed that the object - patient's organs - being scanned is stationary. This assumption is 
violated when there is any motion such as cardiac motion, blood flow, respiratory motion, or 
patient restlessness during the relatively long period of data acquisition, which is between 0.5 
and 2.5 seconds for third generation X-ray CT scanners. The violation of this assumption causes 
motion artifacts, which appear as blurring, doubling and/or streaking in the reconstructed image. 
As is evident, such artifacts substantially impede diagnosis or even lead to erroneous diagnosis. 

[003] Several techniques for removing these artifacts have been disclosed in the prior art. For 
example, Crawford, C.R., King, K.F., Ritchie, C.J., and Godwin, J.D.: ''Respiratory 
Compensation in Projection Imaging using a Magnification and Displacement Model", IEEE 
TMI, 15(3), pp. 327-332, 1996, teach a linear model for the motion, while Ritchie, C.J., Hsieh, 
J., Gard, M.F., Godwin, J.D., Kim., Y, and Crawford, C.R.: "Predictive Respiratory Gating: A 
New Method to Reduce Motion Artifacts in CT Scans", Radiology, 190, pp. 847-852, 1994, 
model the motion as periodic and take projection data at a particular point in the motion cycle, 
producing a stroboscope-like effect. However, the organ motion is much more complex and 
these techniques have a severely limited ability for correcting organ motion. In a more general 
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technique disclosed by Chiu, Y.H., and Yau, S.F.: "Tomographic Reconstruction of Time- 
Varying Object from Linear Time-Sequential Sampled Projections", Proc. IEEE Conf. ASSP, 
Adelaide, 1, pp. V309-V312, 1994, motion effects are iteratively suppressed from the projection 
data. This process requires proper initialization in order to achieve an acceptable convergence 
period. Another technique based on retrospective gating is taught by Morehouse, C.C., Brody, 
W.R., Guthaner, D.F., Breiman, R.S., and Harell, G.S.: "Gated Cardiac Computed Tomography 
with a Motion Phantom", Radiol., 134(1), pp. 213-217, 1980, and employs an ECG signal to 
identify the diastole phase of a patient's cardiac cycle - i.e. when the heart is least active - and 
use the projection data collected during the diastole phase to reconstruct the tomographic image. 

[004] Although the correction methods disclosed in the above cited references reduce the 
adverse effects of motion to some extent, none of these methods substantially removes the 
motion artifacts from the reconstructed tomographic images. 

Summary of the Invention 

[005] It is, therefore, an object of the invention to provide a method and system for 
processing X-ray CT projection data to track organ motion. 

[006] It is further an object of the invention to substantially improve image quality of 
tomographic images by substantially removing motion artifacts. 

[007] It is yet further an object of the invention to provide a method for processing X-ray CT 
projection data for correcting organ motion artifacts for cardiac imaging. 

[008] In accordance with the present invention there is provided a method for tracking motion 
phase of an object comprising: 

receiving a plurality of projection data indicative of at least a cross-sectional image of the object, 
the projection data being in dependence upon a detector signal indicative of attenuation of an X- 
ray beam along a line through the object between an X-ray source and a detector, the projection 
data being acquired such that the projection data are in dependence upon the attenuation along 
different lines through the object and in dependence upon the attenuation along same lines 
through the object at different time instances; 
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processing the projection data for determining motion projection data of the object indicative of 
motion of the object based on the attenuation along at least a same line through the object at 
different time instances; 

selecting a motion phase of the object from the motion projection data; and, 

selecting the projection data acquired at time instances within the selected motion phase of the 

object. 

[009] In accordance with the present invention there is provided a method for tracking motion 
phase of an object comprising the steps for: 

receiving a plurality of projection data indicative of at least a cross-sectional image of the object, 
the projection data being in dependence upon a detector signal indicative of attenuation of an X- 
ray beam along a line through the object between an X-ray source and a detector, the projection 
data being acquired such that the projection data are in dependence upon the attenuation along 
different lines through the object and in dependence upon the attenuation along same lines 
through the object at different time instances; 

processing the projection data for determining motion projection data of the object indicative of 
motion of the object based on the attenuation along at least a same line through the object at 
different time instances; 

selecting a motion phase of the object from the motion projection data; and, 

selecting the projection data acquired at time instances within the selected motion phase of the 

object. 

[0010] In accordance with the present invention there is further provided a method for tracking 
motion phase of an object comprising: 
providing an X-ray source and at least a detector; 

using the at least a detector sensing X-ray beam radiation attenuated by the object along a line 
through the object between the X-ray source and the detector and providing projection data in 
dependence thereupon; 

rotating the X-ray source and at least a detector around the object for acquiring projection data in 
dependence upon the attenuation of the X-ray beam along different lines through the object; 



3 



138-01 US Patent 

rotating the X-ray source and at least a detector around the object a plurality of times for 
acquiring projection data in dependence upon the attenuation of the X-ray beam along same lines 
at different time instances; 

processing the projection data for determining motion projection data of the object indicative of 
motion of the object based on the attenuation along at least a same line through the object at 
different time instances; 

selecting a motion phase of the object from the motion projection data; and, 

selecting the projection data acquired at time instances within the selected motion phase 

of the object. 

[0011] In accordance with the present invention there is further provided a method for tracking 
motion phase of an object comprising the steps for: 
providing an X-ray source and at least a detector; 

using the at least a detector sensing X-ray beam radiation attenuated by the object along a line 
through the object between the X-ray source and the detector and providing projection data in 
dependence thereupon; 

rotating the X-ray source and at least a detector around the object for acquiring projection data in 
dependence upon the attenuation of the X-ray beam along different lines through the object; 
rotating the X-ray source and at least a detector around the object a plurality of times for 
acquiring projection data in dependence upon the attenuation of the X-ray beam along same lines 
at different time instances; 

processing the projection data for determining motion projection data of the object indicative of 
motion of the object based on the attenuation along at least a same line through the object at 
different time instances; 

selecting a motion phase of the object from the motion projection data; and, 

selecting the projection data acquired at time instances within the selected motion phase 

of the object. 

[0012] In accordance with the present invention there is yet further provided a storage medium 
having stored therein executable commands for execution on a processor, the processor when 
executing the commands performing: 
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receiving a plurality of projection data indicative of at least a cross-sectional image of the object, 
the projection data being in dependence upon a detector signal indicative of attenuation of an X- 
ray beam along a line through the object between an X-ray source and a detector, the projection 
data being acquired such that the projection data are in dependence upon the attenuation along 
different lines through the object and in dependence upon the attenuation along same lines 
through the object at different time instances; 

processing the projection data for determining motion projection data of the object indicative of 
motion of the object based on the attenuation along at least a same line through the object at 
different time instances; 

selecting a motion phase of the object from the motion projection data; and, 

selecting the projection data acquired at time instances within the selected motion phase 

of the object. 

[0013] Reconstructed tomographical images clearly show a substantial improvement in image 
quality by successfully removing motion artifacts using the method for tracking motion phase of 
an object according to the invention. Employment of the method for tracking motion phase of an 
object according to the invention is highly beneficial in cardiac imaging using X-ray CT scan. 
The method is easily implemented as retrofit into existing X-ray CT scanners including helical 
X-ray CT scanners, for example, as executable commands provided on a storage medium for 
execution on an existing signal processing system. 

Brief Description of the Figures 

[0014] Exemplary embodiments of the invention will now be described in conjunction with the 
following drawings, in which: 

[0015] Figure 1 is a simplified diagram schematically illustrating parallel beam projection in 
an X-ray CT scanner; 

[0016] Figure 2 is a simplified diagram schematically illustrating fan beam projection in an X- 
ray CT scanner; 

[0017] Figure 3 is a simplified flow diagram illustrating processing of projection data 
employing a method for tracking motion phases of an object according to the invention; 
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[0018] Figures 4a to 4c are diagrams illustrating simulated projection data for a stationary 
object; 

[0019] Figures 5a to 5f are diagrams illustrating the processing of projection data using a first 
embodiment of a method for tracking motion phases of an object according to the invention; 

[0020] Figure 6 is a simplified diagram illustrating geometrically operation of an unwrapping 
filter according to the invention in the first embodiment; 

[0021] Figure 7 is a simplified diagram illustrating conversion of fan beam data into parallel 
beam data; 

[0022] Figure 8 is a simplified diagram illustrating the Shepp and Logan head phantom; 

[0023] Figures 9a to 9d are diagrams illustrating the processing of projection data obtained 
from single detector measurements using the first embodiment of the method for tracking motion 
phases of an object according to the invention; 

[0024] Figures 10a to lOd are diagrams illustrating the processing of projection data obtained 
from sum of detector measurements using the first embodiment of the method for tracking 
motion phases of an object according to the invention; 

[0025] Figures 1 la to 1 Id are reconstructed tomographical images; 

[0026] Figures 12a to 12c are diagrams illustrating the processing of projection data using a 
second embodiment of the method for tracking motion phases of an object according to the 
invention; 

[0027] Figures 13a to 13f are diagrams illustrating the processing of projection data of a first 
patient using the first and the second embodiment of a method for tracking motion phases of an 
object according to the invention; 

[0028] Figures 14a and 14b are reconstructed tomographical images; 
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[0029] Figures 15a to 15e are diagrams illustrating the processing of projection data of a 
second patient using the first and the second embodiment of a method for tracking motion phases 
of an object according to the invention; 

[0030] Figures 16a and 16b are reconstructed tomographical images; 

[0031] Figure 17 is a simplified block diagram illustrating a system implementation of the two 
embodiments of the method for tracking motion phases of an object according to the invention; 

[0032] Figure 18a is a diagram illustrating 360° interpolation of helical projection data; 

[0033] Figure 18b is a diagram illustrating a sinogram corresponding to Fig. 18a; 

[0034] Figure 19a is a diagram illustrating 180° interpolation of helical projection data; 

[0035] Figure 19b is a diagram illustrating a sinogram corresponding to Fig. 19a; 

[0036] Figures 20a and 20b are schematic block diagrams illustrating a single-slice and a 
multi-slice X-ray CT scanner, respectively; 

[0037] Figure 21 is a diagram illustrating a first implementation of the method for tracking 
motion phases of an object with interpolation followed by identification of motion phase; 

[0038] Figure 22 is a diagram illustrating a second implementation of the method for tracking 
motion phases of an object with identification of motion followed by interpolation; 

[0039] Figures 23a to 23d are diagrams illustrating application of the first implementation 
using Constant Attenuation property of Radon Transform CART for identifying motion phase 
after 360° interpolation in comparison with the processing of non-helical projection data; and, 

[0040] Figures 24a to 24d are diagrams illustrating application of the second implementation 
using CART for identifying motion phase before 180° interpolation in comparison with the 
processing of non-helical projection data: 
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Detailed Description of the Invention 

[0041] In X-ray CT, projection data is processed using a tomographic image reconstruction 
process producing a cross-sectional image of a patient's body. For the data acquisition and the 
image reconstruction process it is assumed that the object - patient's organs - being scanned is 
stationary. This assumption is violated when there is any motion such as cardiac motion, blood 
flow, respiratory motion, or patient restlessness during the relatively long period of data 
acquisition, which is between 0.5 and 2.5 seconds for third generation X-ray CT scanners. In 
particular, motion artifacts caused by cardiac motion substantially impede diagnosis or even lead 
to erroneous diagnosis in cardiac imaging. 

[0042] In order to overcome this problem and to enable use of X-ray CT scanners in cardiac 
imaging, a method for tracking motion phase of an object according to the invention has been 
implemented into the processing of the X-ray CT projection data. During a cardiac cycle speed 
of heart motion is not uniform. Within a fraction of the cardiac cycle - diastole phase - the heart 
is less active, i.e. there is less motion, compared to the rest of the cardiac cycle. The method for 
tracking motion phase of an object according to the invention, which will be described in detail 
below, identifies the projection data that have been acquired during the "less active" - diastolic - 
phase of the cardiac cycle. Once these projection data are correctly identified they are fed into 
the tomographic image reconstruction process producing images substantially free of cardiac 
motion artifacts. As will become evident to those skilled in the art, the method for tracking 
motion phase of an object according to the invention is not limited for tracking cardiac motion, 
but is also applicable for tracking phase of other types of motion. 

[0043] In order to provide better understanding of the method for tracking motion phase of an 
object according to the invention, the acquisition of the projection data in X-ray CT scanners and 
the process of tomographic image reconstruction from the projection data will be briefly 
discussed in the following. 

[0044] When an X-ray beam 5 passes through an object 1, as shown in Fig. 1, which is 
modeled as a two dimensional distribution f(x,y) of the X-ray attenuation, total attenuation 
suffered by the X-ray beam is determined as the line integral of the attenuation along the straight 
line 5, which is given by: 



8 



138-01 US 



Patent 



xcos# + >>sin# = t , 
with respect to reference line 3. The total attenuation is then given by: 

00 00 

P(9,t)= \ jf(x,y)5(xcos0 + ysm0-t)dxdy. 

-CO -CO 

Here, S(.„) is the Dirac delta function and P{0,t) is known as Radon Transform of the function 
f(x,y). A projection is formed by combining a set of these line integrals. The simplest 
projection is a collection of parallel line integrals as is given by P(6,t) for a constant 6 and 
varying t - parallel projection. In another type of projection known as fan beam projection, the 
line integrals are determined along fan lines, as shown in Fig. 2, given by: 

x cos(y + a)+ y sin(y + cr) = D sin a . 

The line integrals are then given by: 

00 00 

R(y,o)= J ^ f (x, y)s(x cos(y ^ a) + ysin(y + a)- D sin cr)dxdy. 

-00 -CO 

In third generation CT scanners utilizing fan beam projection, the data acquisition process is 
performed using discrete detectors. These discrete projection measurements R(y n , a m ), 
(n = 1,2,..., N), (m = 1,2,..., M) are the line integrals given by: 

R(r n ><7 m )= l\f(x 9 y)S(xcos{y n +a m )+ysin{y n +a m )- Dsma m )dxdy , 

where the projection angle is y n , and the m -th detector within the detector array is oriented at an 
angle o m with respect to the center of the array. The angular step increment between two 
successive projections is defined by: Ay = In I N , At is the time elapsed between two 
successive projections and TV is the total number of projections acquired during the period 
T = NAt required for one full rotation of the source and detector array around the object f(x 9 y) 
being scanned. 
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[0045] For the tomographic image reconstruction Fourier Slice Theorem or, more frequently, 
filtered back projection is applied. For fan beam projection data a modified version of the filtered 
back projection - known as weighted filtered back projection - is used. 

[0046] The tomographic image reconstruction process of an image f(x,y) from the parallel 
projections P(d,t) is briefly outlined as follows: 

f(x 9 y) = - \Q{0,xcosO + ysmO)de 9 (1) 

^ 0 

where 

Q{O,t)=)s(0,w]M\e J2m "dw, (2) 

-oo 

with S(&, w) being the Fourier transform of the projection P(Oj): 

OO 

s{e,w)= \p{e,ty jlm "dt. 

-00 

Equations (2) and (1) represent the operations of filtering and back projection, respectively, 
hence the name filtered back projection. Details concerning filtered back projection, weighted 
filtered back projection as well as computer implementations of the same are disclosed in Kak, 
A.C. and Slaney, M: "Principles of Computerized Tomographic Imaging", New York, IEEE 
Press, 1999. 

[0047] Referring to Fig. 3, a simplified flow diagram of the method for tracking motion phase 
of an object according to the invention implemented into an X-ray CT data acquisition and image 
reconstruction process is shown. An X-ray source and at least a detector of an X-ray CT scanner 
100 are rotated around the object - box 102. The detector senses the X-ray beam radiation 
attenuated by the object along a line through the object between the X-ray source and the 
detector - box 104 - and provides projection data in dependence thereupon - box 106. During 
one full rotation of the X-ray source and the detector projection data in dependence upon the 
attenuation of the X-ray beam along different lines in a cross-sectional plane through the object 
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are acquired. By rotating the X-ray source and the detector several times around the object, 
projection data in dependence upon the attenuation of the X-ray beam along same lines at 
different time instances are acquired. Using a processor 1 10 the received projection data - box 
1 12 - are processed according to the invention to determine the motion phase of the object with 
the object having the least motion. The projection data in dependence upon the attenuation of 
the X-ray beam along same lines through the object acquired at different time instances are 
compared - box 1 14 - and based on changes of the attenuation - indicating motion of the object 
- the motion phase of the object corresponding to the least changes of the attenuation is 
determined - box 116. The projection data acquired at time instances within the determined 
motion phase of the object are then selected for image reconstruction - box 118. Determination 
of the motion phase and selection of the corresponding data is performed, for example, by a user 
on a graphical display of the processed data, or using a processor for comparing the processed 
data with a threshold. The selected projection data are then processed for tomographic image 
reconstruction using, for example, filtered back projection - box 120 - providing image data 
indicative of a cross-sectional image of the object, which are then displayed - box 122 - using, 
for example a monitor 120. 

[0048] In the following, two embodiments of the method for tracking motion phase of an 
object according to the invention, as indicated above in boxes 1 14 to 1 18, will now be described 
in more detail. In the first embodiment a Spatial Overlap Correlator (SOC) is combined with an 
unwrapping filter for tracking the cardiac motion and accurately detecting the less active phase 
therefrom. The SOC concept has been described in detail by the present inventor in 
Stergiopoulos, S.: "Otpimum Bearing Resolution for a Moving Towed Array and Extension of its 
Physical Aperture ", JASA, 87(5), pp. 2128-2140, 1990, and in Stergiopoulos, S.: 
"Implementation of Adaptive and Synthetic Aperture Processing in Real-Time Sonar Systems", 
Proc IEEE, 86(2), pp. 358-396, 1998. The process of SOC is best understood by a simple 
example. Consider two photographs taken at time instances f 0 and t Q + A^of a field of view 
including a moving train and stationary objects such as trees and houses. If corresponding pixels 
of the two photographs are subtracted from each other, the resulting image includes only 
information about the motion of the train during At . In a similar fashion, projection data acquired 
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at different time instances t 0 and t 0 + At but along same lines through the object are compared to 
detect motion of the object during At . 

[0049] Consider a fixed detector m 0 in the detector array of an X-ray CT scanner. The detector 
is rotated p times around a given stationary object and acquires projection datax[«] . The 
projection data x[n] when plotted versus time n is a periodic signal with period N , which is the 
total number of projections during a full rotation around the object. It is noted that x[n] is zero 
outside the interval 0 < n < pN . A circularly shifted version x s [n] of x[n] is defined as follows: 

r , fxfri + tf] for 0<n<(p-l)N 

XAn *-\ x [n-{p-l)N] for (p-\)N<n<pN { ) 

In other words, the projection data x[n] have been shifted to the left by one period N and 
wrapped around the first period at the end on the right, as shown in Figs. 4a - projection data 
x[n] - and 4b - circularly shifted version x s [n]- for simulated projection data for a stationary 
object with the number of rotations p = 5 and the number of projections per rotation N = 1408 . 
Since x[n]is periodic with period TV it is not possible to detect a difference between x[n] and 
x s \n\. SOC of x[n] is defined as: 

For a given stationary object this signal is zero for all n , as shown in Fig. 4c. For a given 
projection n = n 0 , the two signal values x[n 0 ] and x s [n 0 ] represent the line integral of the given 
object along a same line but at different times. This is because the X-ray source and the detector 
m 0 are disposed at same physical locations at both time instances. If the object is stationary the 
measurements are identical. 

[0050] Now considering the case when a weaker signal y[n] is superimposed on x[n] 
generating a new signal 

z[rc]=jt[«]+j;[4 (4) 
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Fig. 5a illustrates an X-ray CT projection signal x[n] observed by a detector scanning a patient 
with no motion. It is periodic with a period equivalent to a full rotation - 1408 projections per 
second. Fig. 5b is a heartbeat signal y[n] with 95 beats per minute. Fig. 5c illustrates a projection 
signal z[n] observed by the detector when scanning the patient with heartbeat. The goal is to 
recover the signal y[n] from the observed projection signal z[n] - Fig. 5c. In order to recover the 
heartbeat signal y[n] from the observed projection signal z[n], the observed projection signal 
z[n] is first circularly shifted by one period N resulting in z 5 [n], shown is Fig. 5d. Then the 
SOC of the observed signal z[n]is determined as 

z soc M = A n Y z s M = *M+ y[*l- x, W- y* W 

= x[n] - x s [n]+y[n]- y s [n] = y[n]- y s [n] . 

The reason for the term y[n]- y s [n] not to disappear is because it is not periodic with the period 
N . Taking the SOC of the observed signal z[n] cancels the portion of z[n] with period N - 
x[n] - and leaves the portion that is not periodic with period N - the term^[«]- y s [n], shown in 
Fig. 5e. In other words, the stationary signal components have been removed leaving only the 
signal portions due to the moving/deforming object. However, due to the circular shifting of the 
observed signal and wrapping of the first period, information about a time reference is lost in the 
signal shown in Fig. 5e. In the following steps the lost time information is recovered by 
obtaining y[n] from ^ [«]. Considering pN point Fast Fourier Transforms (FFTs) of 

y[n], y s [n], and as follows: 

P N-\ J UlL 

y [ n ]^£L^ Y [k]= 2>[*> [pN) for k = 0X...,pN-\ 

/i=0 

zsoc W =M- ys W* JT] -* z soi W = Y[k 
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Now it is possible to recover ^[/?]from z soc [n] = y[n]- y s [«]as: 



k 



^SOC [fi] 



1-g 



-<7 



1 - cos 



K P 



+ 7 sin 



(6) 



Using the relation in equation (6) the lost information about a time reference in the signal shown 
in Fig. 5e is recovered resulting in the heartbeat signal y[n] shown in Fig. 5f. A comparison with 
the signal shown in Fig. 5b indicates that the heartbeat motion phases have been correctly 
recovered. 

[0051] It is noted that the denominator in equation (6) becomes zero whenever k I p is an 
integer - which will happen N times. In order to avoid "dividing-by-zero", the denominator is 



set to be s + j sin 





( 

when 1 - cos 


2nk^ 






. P J 



< € . It is interesting to note that the denominator of 



equation (6) represents geometrically traversing a circle at discrete points exactly N times, as 
shown in Fig. 6. Further, if the amount of shift to produce z s [n] is one half the total length of 

z[n] - i.e. p = 2 in the current setting - the denominator of equation (6) reduces to 1 - (- 1)* - i.e. 
0,2,0,2,0,.... for k = 0,1,2,3,4,... respectively. This corresponds to the two points indicated by "x" 
in Fig. 6 and simply scales down the FFT of z 50C [n]. Consequently, it is not possible to recover 
y[n] from >{«] in this case. 

[0052] While in the above example z s [n] was generated by shifting z[n] equivalent to one 
rotation - i.e. N projections - it is possible to shift z[n] by half rotation N 12 in case of a middle 
detector of a detector array, because the middle detector measures the line integral along a same 
line after every half rotation - i.e. the signal detected by the middle detector is periodic with 
period N 12. For other detectors of the detector array it is possible to determine half rotation 
SOC by processing the projection data obtained from the detector array - parallel beam sinogram 
data - using the symmetry property of the Radon transform. This results in a significantly lesser 
amount of projection data to be acquired exposing a patient to a lesser amount of radiation. 
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[0053] The processing according to equation (6) is referred to herein as an "unwrapping filter" 
as it restores the temporal information of y[n], which was lost - or wrapped up - in y[n\- y s [n] . 

It is noted that though the temporal contents of y[n] are restored, the DC component of y[n] is 
lost because y[n)- y s [n] is a zero mean signal. In order to account for the effect of shifting and 
wrapping as defined by y s [n] in equation (3) it is important to apply a corresponding FFT - here 
pN point FFT - in equation (6). 

[0054] Summarizing, if a weaker signal y[n] caused by cardiac motion is superimposed on a 
regular projection signal Jt[«]due to a stationary object resulting in an observed projection 
signal z[n] , the SOC processing shown above removes the signal component caused by the 
stationary object. Using the unwrapping filter temporal information lost during the SOC process 
is recovered allowing identification of the less active phase of the cardiac cycle versus the 
temporal axis. Figs. 5c to 5f illustrate the signal processing steps of SOC and unwrapping filter. 

[0055] In the following, a second embodiment of the method for tracking motion phase of an 
object according to the invention will be described. Here, the less active phase of the cardiac 
cycle is identified using a property of the Radon transform. The property will be briefly 
described before its application for tracking the motion phase of an object is explained. 

[0056] Let P(0j)bc the Radon transform of an object f(x,y) contained within a unit circle, 

i 

which is possible without loss of generality. The integral jp(0j)t k dt is a homogeneous 

-i 

polynomial of degree kin cos£ and sin where k = 0,1,... i.e. 



\p(6j)t k dt = Y^aj cos k ~ J 0sin J 0 Vk>0 (k integer). 



7=0 



For k = 1 follows 



i . f 

^P(6j)tdt = a 0 cos0 + a x sin (9 = ^a] + a] sin t 
-i V 



<9 + tan 



-i 



(a, ^ 



\ a \ J 
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From the above equation follows that the center of mass ^P(Oj)tdt has a sinusoidal dependence 

-i 

on 9 . Even in the case of the function f{x,y) consisting of a single point it is mapped to a 
sinusoid in P{0,t) . Therefore, the function P(9,t) is called a "sinogram". 

For k = 0 follows 
i 

\p(9,t)dt = a Q (i.e. constant) V£ (7) 

The constant on the right hand side of equation (7) is called the "mass of the object" and is a 
measure of the total X-ray attenuation constant of the object f(x 9 y) . It is evident that the mass of 
a given object - or the total X-ray attenuation constant - is constant regardless of a projection 
angle from which it is viewed. 

[0057] Equation (7) states that the integral of the Radon transform with respect to t is constant 
for all 9 for a given stationary object. The data acquisition process in parallel beam projection 
results in a 2-D array called "sinogram matrix". The dimensions of this matrix are pN - total 
number of projections - by M - number of detectors in the detector array. According to the 
property of equation (7), the sum of this matrix along M results in a constant vector of length 
pN provided the object was stationary. If the objects starts deforming at some projection, the 
total X-ray attenuation constant starts changing and as a result the vector will start changing at 
the same projection. If the vector starts deviating from its constant value at a projection n Q it is an 
indication of a start of motion of the object at that projection time. Therefore, the motion phase 
of an object is tracked by integrating the Radon transform of the projection data according to 
equation (7) - i.e. the sum along M of the sinogram matrix - and tracking the changes of the sum 
of the sinogram matrix. 

[0058] It is noted that the above equations are valid for parallel beam projection data and not 
for fan beam projection data. In order to apply the method for tracking the motion phase of an 
object as outlined above the fan beam data are re-sorted into equivalent parallel beam data. This 
process is known as "rebinning" and performed according to the following relationship: 
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t = Dsina and 0 = y + cr. 

The parameters of the above relationships are described above and illustrated in Figs. 1 and 2. 
The process of rebinning from fan beam data I to parallel beam data II is shown in Fig. 7. The 
vertical dashed line indicates the location of the central detector. As a result of the rebinning, the 
data corresponding to a single fan projection is spread over a range equivalent to the fan angle in 
the corresponding parallel projection. Due to this spread the effect of a starting motion of the 
object is spread over a range when transformed from fan beam to parallel beam data. 

[0059] In the following, applications of the method for tracking motion phase of an object 
according to the invention will be described using simulated data and real patient data. CT 
scanner projection data for the Shepp and Logan head phantom have been simulated. The 
phantom comprises 10 ellipses, shown in Fig. 8, with one - ellipse "d" - pulsating for simulating 
cardiac motion. Projection data for 5 rotations of the CT scanner have been simulated with 1056 
projections per rotation. The pulsating ellipse was pulsating with 6 pulses during this time 
interval with a period equivalent to 924 projections. The cycle of the pulsating movement has 
been set to 43%, i.e. the ellipse was moving for 43% of the time interval and stationary for 57% 
of the time interval. Implementation of the SOC technique on the simulated phantom data 
presented some problems in precisely identifying the stationary period of the pulsating ellipse, 
shown in Figs. 9a to 9d. Fig. 9a illustrates the simulated projection data of the pulsating phantom 
detected by one detector. Fig. 9b is the corresponding SOC projection signal. Comparison of 
Figs. 9c and 9d, the SOC projection signal after application of the unwrapping filter and the 
motion reference signal of the ellipse "d", respectively, illustrates the difficulty in identifying the 
stationary phase of the phantom in this case. The reason for this inaccuracy is that the line 
corresponding to an X-ray beam between the source and the detector is not intersecting the 
pulsating ellipse during the whole rotation of the detector. In other words, the pulsating ellipse is 
not in the line of view for some source-detector locations during their rotation around the 
phantom object. Hence, equation (5) did not hold during the whole rotation rendering precise 
detection of the stationary period of the pulsating ellipse impossible. This problem is easily 
overcome by ensuring that at any time during the rotation, a source-detector line is intersecting 
the pulsating ellipse by employing, for example, a plurality of detectors instead of a single one 
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and adding the measurements acquired by all the detectors for each projection and forming a 
single signal. 

[0060] Figs. 10a to lOd illustrate the simulated projection signals employing the above 
modified signal processing. In Fig. 10a, the sum of all detector measurements is shown. Figs. 
10b and 10c show the output data from the SOC process and the unwrapping filter, respectively. 
A comparison of the output data of the unwrapping filter with the actual motion phases of the 
pulsating ellipse, shown in Fig. lOd, clearly illustrates that the motion phases are precisely 
identified using the above method for tracking motion phase of an object according to the 
invention. Finally, Fig. 1 la to lid illustrate reconstructed images based on different portions of 
the sinogram data. Fig. 11a shows the image of the stationary phantom for comparison. For the 
current simulation there are 681 projections in 180° + fan angle. Using the projection data 
obtained in the beginning from projection number 0 to 681 in a standard tomographic image 
reconstruction without employing the method for tracking motion phase of an object according 
to the invention will result in an image shown Fig. 1 lb. Here, motion artifacts due to the 
pulsating ellipse are clearly visible. Using the motion phase information - Fig. 10c - a single 
portion of the projection data from projection number 397 to 1078 is selected. As indicated by 
the signal shown in Fig. 10c there is some motion present. Using these not absolutely motion free 
data results in the reconstructed image shown in Fig. 11c. Here, some motion artifacts are present 
but substantially less than in the image shown in Fig. 1 lb. Finally, two portions of the projection 
data are selected comprising motion free data - from projection number 397 to 8 14 and from 
2245 to 2509 - and used for reconstructing the image shown in Fig. lid. Comparison with Fig. 
11a shows substantial removal of the motion effects of the pulsating ellipse. 

[0061] Figs. 12a to 12c illustrate the processing of fan beam projection signals using the 
second embodiment of the method for tracking motion phases of an object according to the 
invention based on the Radon transform property. Fig. 12a shows the sum of detector 
measurements without rebinning. Converting the fan beam projection data into parallel beam 
data using rebinning enhance the motion phase information as shown in Fig. 12b. A comparison 
of the signal shown in Fig. 12b with the actual motion phases of the pulsating ellipse, shown in 
Fig. 12c, clearly illustrates that the motion phases are precisely identified using the second 
embodiment of the method for tracking motion phases of an object according to the invention. 
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Selection of projection data for image reconstruction results in a same image as shown in Fig. 
1 Id. The rebinning process leaves some data points empty in the parallel projection resulting in 
signal distortions at the beginning and the end of the diagram shown in Fig. 12b. The length of 
the distortions is half the fan angle at each end. It is obvious that these distortions do not affect 
the tracking of the motion phases and the projection data selection. 

[0062] In the following, the method for tracking motion phase of an object has been applied to 
projection data obtained in clinical tests. For these tests the X-ray CT scanner was positioned to 
obtain a patient's heart ventricles. The patient was asked to follow breath holding procedures. 
The CT period of rotation was set to 1.0 second - 1408 projections. Projection data were 
recorded for 20 rotations. For this test, the patient's ECG signal was also recorded. The heart 
beat was measured to be 75 beats per minute (bpm). In Fig. 13a 5 the sum of rows of a parallel 
beam sinogram is shown. Half rotation SOC of the signal in Fig. 13a is shown in Fig. 13c. The 
result of the unwrapping filter is depicted in Fig. 13d. Figs. 13a and 13d reveal the motion phase 
information. In order to more clearly identify the motion phase information, the data shown in 
Figs. 13a and 13d are further processed employing a sliding window process step. In Fig. 13b, 
the standard deviation of a sliding window equivalent to 180° + fan angle projections - which are 
908 projections - is shown, i.e. the value shown at a projection n 0 in Fig. 13b is the standard 

deviation of a window from n Q to n 0 + 908 of Fig. 13 a. The same processing has been performed 

with the data shown in Fig. 13d with the results shown in Fig. 13e. A high value at a projection 
point in Figs. 13b and 13e indicates large fluctuations in a data window starting from that point 
predicting larger motion artifacts in the image reconstructed from the projection data within that 
window. Accordingly, a small value indicates less fluctuation in the data window predicting an 
image with lesser motion artifacts. As indicated in Figs. 13b and 13e, the amplitude of overall 
movement increases significantly during the latter part of the CT data acquisition process. 

[0063] Fig. 13f shows the triggering signal based on the patient's ECG, which is used to 
implement ECG retrospective gating for identifying phases of cardiac motion. The square pulse 
of a triggering signal coincides with the QRS peaks of the ECG pulse, which indicates the start 
of diastole phase when the heart is moving the least. Figs. 13b and 13d show a low value at 
projection number 1 1 127 - indicated by the first vertical dashed line - indicating lesser motion 
near the apical portion of the ventricles in a data segment starting at this projection. The image 
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reconstructed from this data segment, shown in Fig. 14a, confirms this, as the heart's 
pericardium and its thickness in the upper part of the image are clearly visible. Also, the circular 
lobe of the liver is visible in the left part of the image. According to Fig. 13f, the data segment 
starting at projection number 1 1 127 includes the QRS phase of the heart cycle suggesting 
strongest heart motion. This is contrary to the observation using the method according to the 
invention and Fig. 14a showing least motion artifacts. Another interesting observation is made at 
projection number 24977 - indicated by the second vertical dashed line. Here, the two 
embodiments of the method according to the invention - Figs. 13b and 13e - indicate large 
motion amplitude whereas the ECG signal suggests less motion in a data segment starting from 
this projection. The image reconstructed from this data segment - Fig. 14b - shows large motion 
artifacts resulting in an image of substantially inferior quality. 

[0064] For a second patient, the CT scanner period was again set to 1 .0 second - 1408 
projections - and projection data were recorded for 28 rotations. During the data acquisition 
process the patient was asked to breathe freely. The patient's heartbeat was recorded to be 72 
bpm. The signal processing steps shown in Figs. 15a to 15e are similar to the ones shown in 
corresponding Figs. 13a to 13e. The sum of rows of a parallel projection beam sinogram is 
shown in Fig. 15a directly providing motion phase information. The output of half rotation SOC 
and unwrapping filter are shown in Figs. 15c and 15d, respectively. Fig. 15d also provides 
motion phase information similar to Fig. 15a. The standard deviation of a sliding window 
equivalent to 180° + fan angle projections applied to the data shown in Figs. 15a and 15d are 
shown in Figs. 15b and 15e, respectively. From the information given in Figs. 15b and 15e, two 
projection points - indicated by vertical dashed lines - have been selected to illustrate the 
capability of the method for tracking motion phases of an object according to the invention. The 
data segment starting from projection number 17429 shows high values in Figs. 15b and 15e 
indicating large motion. The image reconstructed from this data segment shown in Fig. 16a 
confirms the prediction showing large motion artifacts. On the other hand, the data segment 
starting at projection number 16763 shows low values in Fig. 15b and 15e indicating small 
motion. The image reconstructed from this data segment is illustrated in Fig. 16b showing 
substantially less motion artifacts making calcification clearly visible. 
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[0065] The reconstructed images indicate a substantial improvement in image quality by 
successfully removing motion artifacts using the method for tracking motion phase of an object 
according to the invention. 

[0075] Referring to Fig. 17, a system implementation 10 of the two embodiments of the method 
for tracking motion phase of an object according to the invention is shown. For example, 
executable commands for execution on a processor 1 10 of a workstation 14 are stored in non- 
volatile memory 16. When executing the stored commands the processor 110 performs the 
processing steps of one of the above embodiments of the method for tracking motion phases of 
an object. The workstation is linked via port 18 to an X-ray CT scanner 100 for data 
communication therewith. Projection data acquired from an object 12 are transmitted via port 18 
to the workstation 14 for processing. After processing using the method for tracking motion 
phase of an object, selected projection data are used for reconstructing tomographic images, 
which are then, for example, provided to a monitor 120 of the workstation 14 for display to a 
medical practitioner. 

[0076] Furthermore, it is possible to implement the two embodiments of the method for tracking 
motion phases of an object according to the invention as retrofit into existing X-ray CT scanners, 
for example, as executable commands provided on a storage medium for execution on an 
existing signal processing system. 

[0077] According to the invention, the two embodiments of the method for tracking motion 
phase of an object and reconstructing motion artifacts free image are also implementable into the 
data processing of projection data acquired by multi-slice helical CT scanners, as will be 
described later. The processing of the projection data comprises two operations: tracking the 
motion phase and interpolating the missing data. First the operation of interpolation is described. 
In order to simplify the description, only single-slice helical CT scanner measurements are 
considered. As will become evident, the same data processing concept is easily extended to 
process projection data acquired by multi-slice helical CT scanners. The helical CT scan involves 
a simultaneous translational motion of a patient 12 through gantry 101 along an axis of rotation 
of an X-ray source and detectors with the X-ray source and the detectors - not shown in Fig. 1 7 - 
rotating around the patient 12 such that continuous helical projection data are acquired 
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throughout a volume of interest comprising, for example, the patient's heart. The X-ray source 
and the detectors trace a helical path around the patient while detecting projection data. To obtain 
a cross-sectional image, the projection data corresponding to a predetermined cross-section - 
substantially normal to the axis of the translational motion - are determined by interpolating 
projection data acquired along an axis of the translational motion before and after the 
predetermined cross-section. Referring to Fig. 18a, a scan diagram for a 360° interpolation 
process is shown. Slanted lines 200 represent the helical path of the X-ray source and the 
detectors around the patient's body while it is moved through the gantry with respect to 
projection angle 202 and the axis of the translational motion 204. For a predetermined cross- 
section 206 there is only one direct measurement of the projection data available. The other 
projection data corresponding to the cross-section 206 are obtained by interpolating between 
projection data that are 360° apart - or the pitch Az of the helix. These projection data used in the 
interpolation are acquired from a same projection angle but different axial positions before 208 
and after 210 the axial position x of the predetermined cross-section 206. Fig. 18b illustrates the 
same in a corresponding sinogram of projection data obtained by a plurality of detectors 212 
used for the interpolation. In order to generate the projection data of one cross-section 206, i.e. 
projection data corresponding to one 360° rotation, helical projection data acquired during two 
complete rotations - or 720° - of the detector array around the patient's body are used. In the 
case of parallel beam X-ray projection measurements, projection data for only a half rotation - or 
180° - of the detector array is sufficient for image reconstruction. However, in order to 
interpolate the projection data for the cross-section 206 from helical projection data, helical 
projection data acquired during two complete rotations are utilized. In the case of fan beam X- 
ray projection measurements, projection data for (180°+fan) angle are used for image 
reconstruction. Again, helical projection data acquired during two complete rotations are utilized 
to interpolate the projection data for the cross-section 206. 

[0078] In the 360° interpolation process the helical projection data acquired span a range of 720°. 
Interpolation errors are small if the pitch Az of the helix is relatively small. However, for large 
pitch, the interpolation errors become more prominent in the reconstructed image. Employment 
of an 180° interpolation process, shown in Figs. 19a and 19b, reduces these errors. In this 
interpolation process the helical projection data used for the interpolation are only 180° apart. 
According to the symmetry property of the Radon transform, projection data that are acquired 
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180° apart provide the same line integral along a same physical line through an object and hence 
are the same. In the case of helical scanning, helical projection data acquired 180° apart have less 
axial separation compared to helical projection data acquired 360° apart and hence introduce less 
interpolation errors. In order to generate the projection data of one cross-section 206, helical 
projection data acquired during one complete rotation - or 360° - of the detector array around the 
patient's body are used. 

[0079] In the following two implementations of the method for tracking motion phase of an 
object according to the invention for processing helical CT projection data will be disclosed. 
Figs. 20a and 20b illustrate data acquisition using a single slice - Fig. 20a - and a multi-slice - 
Fig. 20b - CT scanner comprising an X-ray source 300 and a single array of detectors 302 or a 
plurality of detector arrays 304, respectively. In modern multi-slice helical CT scanners, four or 
more detector arrays are employed. In order to describe the two implementation methods, it is 
assumed that a distance between two consecutive detector arrays is Az and the speed of the 
translational motion of the patient is such that it covers a distance of Az during one rotation of 
the detector arrays resulting in an overlap of measurements through each cross-section, i.e. each 
helical path is traced by each detector array. 

[0080] In a first implementation of the method for tracking motion phase of an object, the 
interpolation is performed prior motion phase detection. Referring to Fig. 21, acquisition of 
multi-slice helical projection data is shown as a function of translational movement 400 along a 
patient and rotation 402 of four detector arrays 404. For a predetermined cross-section 406 
helical projection data - indicated by thick lines - that are one rotation on each side of the cross- 
section 406 are identified for interpolating the projection data of the cross-section 406. It is noted 
that due to three rotations of the detector arrays 404 and the presence of four detector arrays 404, 
the acquired helical projection data 408 are equivalent to projection data acquired during twelve 
rotations of a single-slice CT scanner. Of the helical projection data corresponding to twelve 
rotations 408, only helical projection data corresponding to six rotations - thick lines - are 
identified for interpolating the projection data of the cross-section 406. The identified helical 
projection data is then used in the 360° interpolation process to generate projection data for the 
cross-section 406. At any given instant, the helical projection data used for interpolating the 
projection data of the cross-section 406 have been acquired by consecutive detector arrays on 
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either side of the predetermined cross-section 406. These helical projection data have been 
acquired at a same time instant but at different locations 400. The interpolated projection data of 
the cross-section 406 is equivalent to projection data generated by a single-slice CT scanner 
scanning for three rotations in one plane, i.e. the interpolated projection data are equivalent to a 
sinogram of three rotations. The interpolated projection data are then processed using one of the 
two embodiments of a method for tracking motion phase of an object to select projection data 
that have been acquired during, for example, the diastole phase of the cardiac cycle. The selected 
projection data are then used for reconstructing a tomographic image corresponding to the cross- 
section 406. 

[0081] In a second implementation of the method for tracking motion phase of an object the 
interpolation is performed after motion phase detection. Referring to Fig. 22, acquisition of 
multi-slice helical projection data is shown as a function of translational movement 500 along a 
patient and rotation 502 of four detector arrays 504. For a predetermined cross-section 506 
helical projection data - indicated by thick lines - half a rotation on each side of the cross-section 
506 are identified for interpolating the projection data of the cross-section 506. It is noted that 
due to four rotations of the detector arrays 504 and the presence of four detector arrays 504, the 
acquired helical projection data 508 is equivalent to projection data acquired during sixteen 
rotations of a single-slice CT scanner. Of the helical projection data corresponding to sixteen 
rotations 508 only helical projection data corresponding to four rotations - indicated by thick 
lines - are identified for interpolating the projection data of the selected cross-section 506. The 
identified helical projection data are then processed using one of the two embodiments of a 
method for tracking motion phase of an object to select projection data that have been acquired 
during a predetermined motion phase, for example, the diastole phase of the cardiac cycle. The 
helical projection data do not precisely belong to the cross-section 506 but have been acquired in 
very close proximity. Therefore, it is possible to treat the helical projection data as belonging to 
the cross-section 506 in the method for tracking a motion phase without introducing substantial 
errors. After selection of the helical projection data belonging to a predetermined motion phase, 
the selected helical projection data are used to interpolate projection data belonging to the cross- 
section 506 using the 180° interpolation process. The interpolated projection data are then used 
for reconstructing a tomographic image corresponding to the cross-section 506. 
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[0082] In the first implementation - 360° interpolation - the multi-slice helical projection data 
used in the interpolation process have been acquired during two rotations of the detector arrays 
compared to one rotation in the second implementation - 180° interpolation. As a result, the 
interpolation errors are larger in the first implementation using 360° interpolation. Employing the 
360° interpolation, the multi-slice helical projection data used in the interpolation process have 
been acquired at a same time instance, whereas in the 180° interpolation the multi-slice helical 
projection data used in the interpolation process have been acquired at different time instances, 
resulting in some motion artifacts in the reconstructed image for a non-stationary object. Using 
180° interpolation helical projection data corresponding to four rotations are available for motion 
phase identification compared to three rotations when using 360° interpolation. Obviously, a 
larger amount of projection data has a larger probability of containing a predetermined motion 
phase. 

[0083] Alternatively, 360° interpolation is used in the first implementation and 180° interpolation 
in the second implementation. 

[0084] In the following the two implementations outlined above will be compared using 
simulated projection data for a 3D phantom containing a plurality of ellipsoids. First, non-helical 
projection data have been generated for a stationary phantom. Employment of the method for 
tracking motion phase of an object on the simulated projection data for the stationary phantom 
shows no artifacts, i.e. a non-varying signal, illustrated in Figs. 23a and 24a. Next, non-helical 
projection data have been generated for the phantom with one ellipsoid pulsating. Employment 
of the method for tracking motion phase of an object correctly identifies the pulsating and non- 
pulsating phases of the phantom, shown in Figs. 23b and 24b. In a third experiment, helical 
projection data have been generated for a stationary phantom. Fig. 23c shows a signal obtained 
after employment of the method for tracking motion phase of an object using the first 
implementation - 360° interpolation followed by motion phase identification. Though the 
phantom was stationary, the signal after motion identification shows some periodic disturbance 
with a period equivalent to one rotation, as illustrated in Fig. 23c. A periodic disturbance of the 
signal - instead of a straight horizontal line as in Fig. 23a - indicates the presence of artifacts due 
interpolation errors. Fig. 24c illustrates a signal obtained after motion phase identification 
applied to the identified helical projection data prior the step of interpolation in the second 



25 



138-01 US 



Patent 



implementation. Again, the signal contains a periodic disturbance indicating the presence of 
artifacts due to non-planar acquisition of the projection data. A comparison shows that the 
amplitude of the periodic disturbance in the signal shown in Fig. 24c - second implementation - 
is smaller than the disturbance in the signal shown in Fig. 23c - first implementation. In a fourth 
experiment, helical projection data have been generated for the phantom having one ellipsoid 
pulsating. Fig. 23 d illustrates a signal obtained after motion phase identification using the first 
implementation. The pulsating and non-pulsating phases have been clearly identified. At the 
same time artifacts due to interpolation errors are visible. A tomographical image reconstructed 
using projection data corresponding to the interval indicated by an arrow in Fig. 23d will have 
motion artifacts substantially removed and contain only minor artifacts due to interpolation 
errors. Fig. 24d illustrates a signal obtained after motion phase identification applied to the 
identified helical projection data prior interpolation in the second implementation. Again, the 
pulsating and non-pulsating phases have been clearly identified but artifacts due to non-planar 
acquisition of the projection data are also visible. Since more projection data are needed for the 
180° interpolation in the second implementation than for the 360° interpolation in the first 
implementation, it was not possible to obtain an interval of contiguous motion free projection 
data within the duration of a diastole phase - indicated by the arrow in Fig. 24d. As a result, a 
reconstructed tomographical image will contain some motion artifacts. Artifacts due to non- 
planar acquisition of the projection data are expected to be insignificant compared to motion 
artifacts of a non motion corrected image. 

[0085] Numerous other embodiments of ,the invention will be apparent to persons skilled in the 
art without departing from the spirit and scope of the invention as defined in the appended 
claims. 
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